Spectral solutions to stochastic models of gene expression with bursts and regulation 
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Signal-processing molecules inside cells are often present at low copy number, which necessitates 
probabilistic models to account for intrinsic noise. Probability distributions have traditionally been 
found using simulation-based approaches which then require estimating the distributions from many 
samples. Here we present in detail an alternative method for directly calculating a probability 
distribution by expanding in the natural eigenfunctions of the governing equation, which is linear. 
We apply the resulting spectral method to three general models of stochastic gene expression: a single 
gene with multiple expression states (often used as a model of bursting in the limit of two states), 
a gene regulatory cascade, and a combined model of bursting and regulation. In all cases we find 
either analytic results or numerical prescriptions that greatly outperform simulations in efficiency 
and accuracy. In the last case, we show that bimodal response in the limit of slow switching is not 
only possible but optimal in terms of information transmission. 



Signals are processed in cells using networks of in- 
teracting components, including proteins, mRNAs, and 
small signaling molecules. These components are usually 
present in low numbers [TJ [S] [3J HI |5J [5] , which means the 
size of the fluctuations in their copy counts is comparable 
to the copy counts themselves. Noise in gene networks 
has been shown to propagate [7] , and therefore explicitly 
accounting for the stochastic nature of gene expression 
appears important when predicting the properties of real 
biological networks. 

Although summary statistics such as mean and vari- 
ance are sometimes sufficient for answering questions of 
biological interest (8] , calculating certain quantities, such 
as information transmission J5J [TO] ED HI EH E3] or es- 
cape properties [THl ESI EH ESI EH], requires knowing 
the full probability distribution. Full knowledge of the 
probability distribution can also be used to discern dif- 
ferent molecular models of the noise sources based on 
recent exact measurements of probability distributions 

OEDEaEI]. 

Much analytical and purely computational effort has 
gone into detailed models of noise in small genetic 
switches [8] EH [25| [26]. The most general description 
is based on the master equation describing the time evo- 
lution of the joint probability distribution over all copy 
counts |27j . Some progress has been made by applying 
approximations to the master equation. For example, a 
wide class of approximations focuses on limits of large 
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concentrations or small switches [SI [2H [2S] . Approxima- 
tions based on timescale separation of the steps of small 
signaling cascades have also been successfully used to cal- 
culate escape properties [TH] HH] [30] . More often, model- 
ers resort to stochastic simulation techniques, the most 
common of which is based on the varying-step Monte 
Carlo method [3TJ [32] ■ Probabilistic modeling of stochas- 
tic systems by simulation requires a computational chal- 
lenge (generating many sample trajectories) followed by 
an even more difficult statistical challenge (parameteriz- 
ing or otherwise estimating the probability distribution 
from which the samples are drawn) [33 . 

In a recent paper |34j we introduce a new method for 
calculating the steady state distributions of chemical re- 
actants, which we call the spectral method. The pro- 
cedure relies on exploiting the natural basis of a sim- 
pler problem from the same class. The full problem is 
then solved numerically as an expansion in the natu- 
ral basis [IS]. In the spectral method we use the ana- 
lytical guidance of a simple birth-death problem to re- 
duce the master equation for a cascade to a set of alge- 
braic equations. We break the problem into two parts: 
a parameter-independent preprocessing step, and the 
parameter-dependent step of obtaining the actual prob- 
ability distributions. The spectral method allows huge 
computational gains with respect to simulations. In prior 
work [33] we illustrate the method in the example of gene 
regulatory cascades. We combine the spectral method 
with a Markov approximation, which exploits the obser- 
vation that the behavior of a given species should depend 
only weakly on distant nodes given the proximal nodes. 

In this paper we expand upon the application of the 
spectral method to more biologically realistic models of 
regulation: (i) a model of bursting in gene expression 
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and (ii) a model that includes both bursts and explicit 
regulation by binding of transcription factor proteins. In 
both cases we demonstrate how the spectral method gives 
either analytic results or reduced algebraic expressions 
that can be solved numerically in orders of magnitude 
less time than stochastic simulations. 

We begin with a model of a multi-state birth-death 
process, a special case of which has been used to describe 
transcriptional bursting [2"2l l3"S"] . We illustrate how the 
spectral method reduces the model to a simple iterative 
algebraic equation, and in the appropriate limiting case 
recovers the known analytic results. We also use this sec- 
tion to introduce the basic notation used throughout the 
paper. Next we explore the problem of gene regulation in 
detail. The main idea behind the spectral method is the 
exploitation of an underlying natural basis for a problem 
which we can solve exactly. We explore four different 
spectral representations of the regulation model used in 
previous work 34J that arise from four natural choices 
of eigenbasis in which to expand the solution (cf. Sec. 
II A and Fig. |2| . All representations reduce the master 
equation to a set of linear algebraic equations, and one 
admits an analytic solution by virtue of the tridiagonal 
matrix algorithm. We compare the efficiencies of the rep- 
resentations' numerical implementations and show that 
all outperform simulation. Lastly, we apply the spec- 
tral method to a model that combines bursting and reg- 
ulation. We obtain a linear algebraic expression that 
permits large speedup over simulation and thus admits 
optimization of information transmission. Optimization 
reveals two types of solutions: a unimodal response when 
the rates of switching between expression states are com- 
parable to degradation rates, and a bimodal response 
when switching rates are much slower than degradation 
rates. 



I. BURSTS OF GENE EXPRESSION 

We first consider a model of gene expression in which 
a gene exists in one of Z stochastic "states," i.e. pro- 
tein production obeys a simple birth-death process, but 
with a state-dependent birth rate. In the special case 
of Z = 2, this corresponds to a gene existing in an on- 
or an off-state due, for example, to the binding and un- 
binding of the RNA polymerase. Such a model has been 
used to describe transcriptional bursting [22 , 35 , and we 



specialize to this case in Sec. I C 



For a general Z-state process, the master equation 
9zP z n -i + (n+ l)K+i - (9z + »K + X) n »'*£' W 



p;< 



describes the time evolution of the joint probability dis- 
tribution p*, where z specifies the state (1 < z < Z), n is 
the number of proteins, g z is the production rate in state 
z, £l zz ' is a stochastic matrix of transition rates between 
states, and jf n denotes differentiation of the probability 
distribution with respect to time. Time and all rates 



have been nondimensionalized by the protein degrada- 
tion rate. Note that conservation of probability requires 



(2) 



The relationship between the transition rates in il zz i 
and the probabilities tt z = ^2 n Pn of being in the zth 
state can be seen by summing Eqn. [T] over n; at steady 
state one obtains 



Q zz '7iv = 0, 

z' 

and normalization requires 



(3) 



(4) 



In the following sections, we introduce the spectral 
method and demonstrate how it can be used to solve 
for the full joint distribution jf n . 



A. Notation and definitions 

We begin our solution of Eqn. [T] by defining the gener- 
ating function |27j G z (x) — ^2 n Pn xU over complex vari- 
able x [49] (note that superscript z is an index, while 
superscript n on x n is a power). It will prove more 
convenient to rewrite the generating function in a more 
abstract representation using states indexed by protein 
number In), 



|G Z > = $»>, 



with inverse transform 

Pn = (n\G z ). 



(5) 



(6) 



The more familiar form can be recovered by projecting 
the position space (x| onto Eqn. [5| with the provision 
that 



(x\n) = x n . 
Concurrent choices of conjugate state 



1 



(n\x) = — — 



and inner product 



dx 



(/I/') = f ^{f\x){x\f) 
2m 



(7) 



(8) 



(9) 



ensure orthonormality of the states, (n\n'} — 6 nn i, as can 
be verified using Cauchy's theorem, 
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where the convention 0(0) = is used for the Heaviside 
function. 

With these definitions, summing Eqn.[T]over n against 
|n) yields 



\G Z ) = -(a+ ~ l)(a" - g z )\G z ) + ^ n zz ,\G z 



(11) 



where the operators a + and a raise and lower protein 
number, respectively, i.e. 



a + \n) 



a \n) 

with adjoint operations 



|n + l>, 
n\n — 1) 



(n\a = (n — 1|, 
(n\a~ = (n + 1|. 



(12) 
(13) 



(14) 
(15) 



As in the operator treatment of the simple harmonic os- 
cillator in quantum mechanics, the raising and lower op- 
erators satisfy the commutation relation [a~,a + ] = 1, 
and a + a~ is a number operator, i.e. a + a~\n) = n\n) |50j . 
This operator formalism for the generating function was 
introduced in the context of diffusion independently by 
Doi j3S] and Zeldovich [37] , and developed by Peliti [3H] . 
A review by Mattis and Glasser introduces and dis- 
cusses the applications of the formalism for diffusion. 

The factorized form of the birth-death operator in Eqn. 
[TT] suggests the definition of shifted raising and lowering 
operators 



b+ 
ir 



making Eqn. 11 



\G X ) = -b + b7\G x 



^zz'\G z 



(16) 
(17) 



(18) 



Since b + b z is a new number operator, it is clear that 
the eigenvalues of the birth-death operator —b + b~ are 
nonpositive integers, i.e. 



b + b7\j z 



J\Jz 



(19) 



where nonnegative integers j index (z-dependent) eigen- 
functions \j z ). In position space the jth eigcnfunction 
is 



(x\j z ) = (x-iye^ x - 1 \ 



with conjugate 



e -g z (x-l) 



(20) 



(21) 



such that orthonormality (j z \j' z ) — Sjj> is satisfied under 
the inner product in Eqn. [9] Note that b+ and b~ raise 



and lower eigenstates \j z ) as in Eqns. I2p5 i.e. 

b + \3z) = \(J + l)z), (22) 

k\jz) = j\U-l) z ), (23) 

(Jz\b + = ((j-lU (24) 

(Jz\k = (j + l)((j + l) z \. (25) 

As we will see in this and subsequent sections, going 
between the protein number basis |n) and the eigenbasis 
\j z ) requires the mixed product (n\j z ) or its conjugate 
(jz\n). There are several ways of computing these ob- 
jects, as described in Appendix A. Notable special cases 
are 



((i = o)» = 1, 

(n\(j = 0) z ) 



(26) 
(27) 



where the latter is the Poisson distribution. 



B. The spectral method 

We now demonstrate how the spectral method exploits 
the eigenfunctions \j z ) to decompose and simplify the 
equation of motion. Expanding the generating function 
in the eigenbasis, 



\Gz)=J2 G > 



\Jz 



(28) 



and projecting the conjugate state (j z \ onto Eqn. 18 
yields the equation of motion 

G) = - 3 G) + n.* E G r & \& ) ( 29 ) 



for the expansion coefficients G| (where the dummy in- 
has been changed to f in Eqn. 29 1 



dex j in Eqn. 
Using Eqns. [9] 
uates to 



20 



and 21 the product (j z \f z >) eval- 



3-3 



e(j-j' + i), 



(30) 



where A 2Z < = g z — g 7 j. Noting that (j z \jz') = 1 and that 
(jz\j' z ) = for j' < j, Eqn. [29] becomes 



G z ; 



-E^E^^ "' 



j'<3 



(j-j'Y- 



(31) 



The last term in Eqn. 31 makes clear that each jth term 
is slaved to terms with j' < j, allowing the to be 
computed iteratively in j. The lower-triangular struc- 
ture of the equation is a consequence of rotating to the 
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eigenspace of the birth-death operator; this structure was 
not present in the original master equation. 
At steady-state, G* obeys 



j'<3 



(j-rv- 



(32) 



from which it is clear that G| can be computed succes- 
sively in j. Since 



G z = {(j = 0) z \G z ) = '£p* n {(j = 0)t 



(33) 



(cf. Eqn. 26), the computation is initialized using Eqns. 
[31141 i.e. 



^n„,Gt = o, 

z' 

E G o = I- 



(34) 
(35) 



Recalling Eqns. [6] and [28] the probability distribution is 
retrieved via 



P 



« = E G >i^> 



(36) 



where the mixed product (n\j z ) can be computed as de- 
scribed in Appendix A. 

There is an alternative way to decompose the master 
equation spectrally. Instead of expanding the generating 
function in eigenfunctions \j z ), which depend on the pro- 
duction rates g z in each state, we may expand in eigen- 
functions parameterized by a single rate g, i.e. 



icy = E G fi-?>< 



where 



<*li> 



(x — l) J e 

e -f(x-l) 
(x - + ! ' 



(37) 

(38) 
(39) 



The parameter g is arbitrary and thus acts as a "gauge" 
freedom. 

We may now partition the birth-death operator as 



- b + K = -b + b- - b + T z 



(40) 



where b = a — g such that the |j) are the eigenstates 
of b + b~ , i.e. 



b + b~\j)=j\j) 



(41) 



and T z = g — g z describes the deviation of each state's 
production rate from the constant g. 



Projecting the conjugate state (j\ onto Eqn. 18 and 
using Eqn. [37] (with dummy index j changed to j ) gives 

G ) = -3G) + E0'|6 + r,|j'>G£ + E . (42) 

j> z' 



Recalling Eqn. [24] Eqn. [42] at steady state becomes 



E 



{n zzl -j5 zz ,)G) 



T Z G 



3-1- 



(43) 



Eqn. [43] is subdiagonal in j, meaning computation of the 
jth term requires only the previous (j — l)th term and 
the inversion of the Z-by-Z matrix (i7 — jl) (where I 



is the identity matrix). It is initialized with Eqns. 34 35 
and solved successively in j . The probability distribution 
is retrieved via 



Pn 



E G l 

3 



n\j), 



(44) 



where (n\j) is computed as described in Appendix A. 

As an example of a simple computation employing the 
spectral method, Fig. [l] shows probability distributions 
for the case of Z = 3 states, corresponding to a gene 
that is either off, producing proteins at a low rate, or 
producing proteins at a high rate. For simplicity we set 
the rates of switching among all states equal to a constant 
lu, making the stochastic matrix 



ft 




(45) 



As seen in Fig. [T] when w« 1 (corresponding to a switch- 
ing rate much slower than the degradation rate) the dwell 
time in each expression state lengthens. The slow switch- 
ing gives the protein copy number time to equilibrate in 
any of the three expression states, resulting in a trimodal 
marginal distribution p n . When w>l (corresponding to 
a switching rate much faster than the degradation rate), 
the system switches frequently among the three expres- 
sion states, resulting in an average production rate. In 
this limit, the expression state equilibrates on a faster 
timescale than the protein number state. 



C. The on/off gene 

For the special case of Z = 2 states, as when a gene 
is either "on" (z — +) or "off" (z = — ), it is useful to 
demonstrate how the spectral method reproduces known 
analytic results [22] [35] . The probability distribution can 
be written in vector form as 



Pn 



(46) 



and defining w + and cj_ as the transition rates to and 
from the on-state, respectively, the stochastic matrix 
takes the form 



n = 

Note that Eqn. [3] implies 



7T_ 
7T-1 



0J_ 
UJ+ 



(47) 



(48) 
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where 
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FIG. 1: An example of the spectral method for Z — 3 states. 
Dotted curves show the joint distribution for each of z = 
1, 2, and 3, and solid curves show the marginal distribution 
p„. The stochastic transition matrix is given in Eqn. |45| and 
the setting of u in each panel is indicated in the upper-right 
corner. Production rates are g\ = 0, gi = 3, and g% = 12 
for all panels. Distributions are calculated via the spectral 
decomposition in Eqn. |43| with g — (g z ) = 5. 



which makes clear that increasing the rate of transition 
to either state increases the probability of being in that 
state. 

From Eqn. [32] the spectral expansion coefficients obey 



^G±-o, ± Gj=o, ± £G^^_; (49 ) 
j'<j u J '' 



where A = A^ = — A_ + . Initializing with Gq = 

u)±/(uj+ + u)—) and computing the first few terms reveals 
the pattern 



mri (i'+^) 



r- n-4 (j" + ^+ + ^- + i) 



(50) 



(51) 



where in the second line the products are written in terms 
of the Gamma function. 

For comparison with known results 22, 35 we write 
the total generating function |G) = J2± \G±) in position 
space, 

G(x) = E<*l G ±>=EE<^±>0±l G ±> (52) 

± ± 3 

= EEo*- 1 )^*" 1 ^ ( 53 ) 



± 3 

U± p fl±(x-l) 



E ' 



$[a,/3;y] = E 



T(j + a) T(J3) yi_ 
r(a) r(j + /3)i! 



(55) 



is the confluent hypergeometric function. As shown in 
Appendix B, in the limit (/_ =0, Eqn. |54| reduces to 



G(x) = + uo-\g+{x - 1)], 

and the marginal p n is given by 



(56) 



n\ r(w + ) r(n + w + +w_) 
x$[w + + n, w + + w_ + n; -,g + ], 



(57) 



in agreement with the results of Iyer-Biswas et al. |35j and 
Raj et al. [25]. We remind the reader that in addition to 
reducing to known results in the special case of Z = 2 
states with a vanishing off-rate, the spectral method is 
valid for any number of states with arbitrary production 
rates. 



II. GENE REGULATION 

Next we consider a two-gene regulatory cascade, in 
which the production rate of the second gene is a func- 
tion of the number of proteins of the first gene. As shown 
in previous work |34j , a cascade of any length can be re- 
duced to such a generalized two-dimensional system using 
the Markov approximation, which asserts that the proba- 
bility distribution for a given node of the cascade should 
depend only weakly on the probability distributions of 
the distant modes given the proximal nodes. 

In the present section, we consider only the general- 
ized two-dimensional equation and explore different ap- 
proaches to solving it. The equation describes two genes, 
each with one expression state, with regulation encoded 
by a functional dependence of the downstream protein 
production rate on the upstream protein copy number. 
In section |IID| we make an explicit connection between 
the on/off gene discussed in section I C and the case when 
the functional dependence is a threshold. Finally, in sec- 
tion [ITT] we combine the two types of models and consider 
a system with regulation and bursts. 



c*[to T ,w + +w_ + l;TA(aj-l)] > (54) 



A. Representations of the master equation 

1 . The \n, m) basts 

Defining n and m as the numbers of proteins produced 
by the first and second gene, respectively, the master 
equation describing the time evolution of the joint prob- 
ability distribution p nm is |34] 

Prim = g n -lPn-l,m + ( n + l)Pn+l,m ~ (Sb + n )Pnm 
+P [q n Pn,m-l + (m + l)p n , m +l - (<?« + m)p nm ] . (58) 
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-II*)- 
[9'«1] 



|j' k j) 



FIG. 2: Summary of the bases discussed in Sec. II A (black) 
and their gauge freedoms (barred parameters in grey) . In the 
top row, neither the m nor the n sector is expanded in an 
eigenbasis; in the middle row, one sector is expanded; and in 
the bottom row, both sectors are expanded. The \j, k) basis 
can be viewed as a special case of the \j,kj) or \j,k n ) basis 
with <jj = q or q n = q respectively. The \j, m) basis is not 
discussed as it is not useful in simplifying the problem. 



generating function [27] G(x,y) = J2 nm Pn m x n y m over 
complex variables x and y, or, in state notation, 



p nm \n,m) 



with inverse transform 



(n, m\G). 



(61) 



(62) 



Summing Eqn. 58 over n and m against \n,m) and em- 
ploying the same operator notation as in Eqns. [T2|T5" 
yields 



where 



and 



\G) = -H\G), 



H = b+b n (n) + pb+b m (n), 



(63) 



(64) 



The function q n describes the regulation of the second 
species by the first, and the function g n describes the 
effective autoregulation of the first species, due either to 
a non-Poissonian input distribution or to effects further 
upstream in the case of a longer cascade [34]. Time is 
rescaled by the first gene's degradation rate, so that each 
gene's production rate (g n or q n ) is normalized by its 
respective degradation rate, and p is ratio of the second 
gene's degradation rate to that of the first. We impose no 
constraints on the form of g n or q n — they can be arbitrary 
nonlinear functions. 



Summing Eqn. 58 over m gives a simple recursion re- 
lation between g n and p n at steady state, from which 
explicit relations are easily identified. If p n is known, g n 
is found as 



(n+1) 



Pn+1 
Pn 



If on the other hand g n is known, p n is found as 



Pn 



Po 



n< 

n'=0 



(59) 



(60) 



with pq set by normalization. Note that if the first species 
obeys a simple birth-death process, g n = g = constant, 



and Eqn. 60 reduces to the Poisson distribution. 

In the current representation (Eqn. 58 1, which we de- 
note the |n, m) basis, finding the steady state solution for 
the joint distribution p nm means finding the null space 
of an infinite (or, effectively for numerical purposes, very 
large) locally banded tridiagonal matrix. More precisely, 
defining N as the numerical cutoff in protein number n or 
m, the problem amounts to inverting an 7V 2 -by-JV 2 ma- 
trix, which is computationally taxing even for moderate 
cutoffs N. 

In order to solve Eqn[58]more efficiently we will employ 
the spectral method. We begin as before by defining the 



b n = «n 


-1, 


(65) 


S+ - a+ 


- 1, 


(66) 


b-(n) = a~ 


- 9n, 


(67) 




- q n - 


(68) 



Here the regulation functions have been promoted to op- 
erators obeying g n \n) = g n \n) and q n \n) — q n \n), sub- 
scripts on operators denote the sector (n or m) on which 
they operate, and the arguments of b~ and 6~ remind us 
that both are n-dependent. 

Eqn. [64] makes clear that if not for the n-dependence 
of the operators the Hamiltonian H would be diago- 
nalizable, or, equivalently, if g n and q n were constants 
the master equation would factorize into two individual 
birth-death processes. We may still, however, partition 
the full Hamiltonian as 



H 



Ho + -Hi, 



(69) 



where Ho is a diagonalizable part (and Hi is the corre- 
sponding deviation from the diagonal form) , and expand 
\G) in the eigenbasis of Ho to exploit the diagonality. 

As with the multi-state system in Sec. |IB[ where we ex- 
pand the solution in two different bases, there are several 
natural choices of eigenbasis of Ho- Fig. [2] summarizes 
these choices diagrammatically: starting in the \n,m) 
basis (at the top of Fig. |2j, one may expand in eigen- 
functions either the first species, yielding the \ j,m) basis 
(left), or the second species, yielding the \n,k n ) basis 
(right; in general we allow the parameter defining the 
second species' eigenfunctions to depend on n to reflect 
the regulation of the second species by the first From ei- 
ther the \j,m) or the \n,k n ) basis, one may expand in 
eigenfunctions the remaining species, yielding either the 
\j,k n ) basis (bottom left), in which the second species' 
eigenfunctions depend on the first species' copy number 
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n, or the \j, kj) basis (bottom right), in which the sec- 
ond species' eigenfunctions depend on the first species' 
cigcnmode number j. Both bases reduce to the \j,k) ba- 
sis (bottom center) when the parameter of the second 
species' eigenfunctions is a constant. 

The \n,m) and \j,m) bases are less numerically useful 
than the other bases: as discussed above and detailed in 
Fig. m the | n, m) basis is numerically inefficient; and the 
| j, m) basis does not exploit the natural structure of the 
problem, since, unlike the other bases, it neither retains 
the tridiagonal structure in n nor gains a lower trian- 
gular structure in k (see the sections below). Each of 
the remaining bases, however, has preferable properties 
in terms of numerical stability and ability to represent 
the function sparsely yet accurately (either using a few 
values of n in the \n, k n ) basis or a few values of j in the 
\j, k), \j, kj), or \ j, k n ) bases, for example). Moreover, the 
equation of motion simplifies differently in each of the dif- 
ferent bases. We present the derivations of the equations 
of motion in the following sections, beginning with the 
\j, k) basis, generalizing to the \j, kj) basis, moving to 
the |n, k n ) basis, and ending with the \j, k n ) basis. 



2. The \ j, k) basis 

For expository purposes we start by recalling the spec- 
tral representation used in previous work [34] . We choose 
the diagonal part of the Hamiltonian to correspond to 
two birth-death process with constant production rates 
9 and q, 



H Q = b nK + P b nX 



(70) 



where b~ = a~ — g and 6~ = a~ — q. As in Sec. IB 
the parameters g and q act as gauge choices: their values 
can be set arbitrarily and thus can affect the numerical 
stability of the method. The nondiagonal part, 



Hi 



(71) 



captures the deviations T n — g — g n and A„ = q — q n 
of the regulation functions from the constant rates. We 
expand the generating function as 



\G) 



jk 



Gjk\j, k), 



(72) 



where \j, k) is the eigenbasis of H , i.e. 

H \j,k) = (j + pk)\j,k). (73) 

The eigenbasis is parameterized by the rates g and q, 
meaning in position space (x\j) is as in Eqn. 38 and sim- 
ilarly for (y\k) with x — > y, j — > k, and g — > q. 

With Eqns. 70 72 projecting the conjugate state (j, k\ 
onto Eqn. [63] yields 

Gjk = Hj + pk)G ]k -J2(j\b^n\j')(k\k')G rkl 

j'k' 

-pY,(k\b+\k')m n \j')G rk , (74) 

j'k' 



(where the dummy indices j and k in Eqn. 72 have been 
changed to j' and kl respectively in Eqn. 74). Recalling 



Eqn. [24] and restricting attention to steady state, Eqn. 
l74l becomes 



= -{j + pk)G jk - T j-i,j'Gfk 
j' 

j' 



(75) 



where the deviations have been rotated into the eigenba- 
sis as 

= (#„|/> =5>1n)(s- 5 „)(n|/), (76) 

n 

Ajj, = {j\A n \j')=J2^\ n W-qn){n\j'}. (77) 

n 

Eq n. |75| is subdiagonal in k and is therefore similar to 
Eqn. |43| in that the last term acts as a source term. It is 
initialized using 

Gj = {j,k = 0\G) = J2Pnm(j\n)(0\m) (78) 

i vm 



(cf. Eqn. 26 1 with known p n (cf. Eqn. 60 1, then solved at 

can 



each subsequent k using the result for k — 1. Eqn. 75 
be written in linear algebraic notation as 



G k = -p{D k + S-T) AGfc-i, (80) 

where G k is a vector over j, bold denotes matrices, and 
Djj, = (j + pk)8jji and Sjj, = 5j-xj> are diagonal and 
subdiagonal matrices, respectively. Eqn. [80] makes clear 
that the solution involves only matrix multiplication and 
the inversion of a J-by- J matrix K times, where J and 
K are cutoffs in the eigenmode numbers j and fc, re- 
spectively. In fact, if the first species obeys a simple 
birth-death process, i.e. g n = g = constant, setting g — g 
makes Tjj' = 0, and, since D fc is diagonal, the solution 
involves only matrix multiplication. The decomposition 
of the master equation into a linear algebraic equation 
results in huge gains in efficiency over direct solution in 
the \n, m) basis; the efficiency of all bases presented in 
this section is described in Sec. IIIBI and illustrated in 
Fig. [3] 



Recalling Eqns. 62 and 72 the joint distribution is re- 
trieved from Gj k via the inverse transform 



Pnm — 



jk 



n\j}G jk (m\k), 



(81) 



a computation again involving only matrix multiplica- 
tion. The mixed product matrices (n\j) and (m\k) are 
computed as described in Appendix A. 
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3. The \j, kj) basis 

The \j, fc) basis treats both genes similarly by expand- 
ing each around a constant production rate. We may in- 
stead imagine an eigenbasis that more closely reflects the 
underlying asymmetry imposed by the regulation, and 
make the basis of the second gene a function of that of 
the first. That is, we expand the first gene in a basis \j) 
with gauge g as before, but now we expand the second 
gene in a basis \kj) with a j-dependent local gauge q~j. 
We write the generating function as 



jk 



and Eqn. 63 at steady state becomes 

jk 







\j,kj) 



(82) 



(83) 



where we have partitioned the Hamiltonian for each j as 



Ho(j) = b+b-+pb+b-{j), 
H x {j) = b+T n + pb+A n (j), 



(84) 
(85) 



with b n = a n — g and f n = g — g n as before and now 

Knti) = Kn ~ ?j and ^n(j) = ?j ~ Qn- Note tnat tnis 

basis enjoys the eigenvalue equation 



H (j)\j,kj) = (j + pk)\j,kj). 



(86) 



Projecting the conjugate state (j, kj\ onto Eqn. 83 
yields, after some simplification (cf. Appendix C), the 
equation of motion 



={j + pk)G jk + Y,Tj-i,3>G. 



j'k 



t=i j> 

where Tjji is as in Eqn. [76] and 

A, r = ( J \A n (j)\f)=Y / (j\n)(q J -qn)(n\f), 

n 

{-Qjj'Y 



(87) 



* 33' 



(89) 



with Qjji — qj — qj/. Eqn. 87 can be written linear alge- 
braically as 



G k = - CD" 



S'T) 



^[(S-^^V^ + pA^V^ 1 ] G fe _ 



(90) 



where Djj, and S^, are defined as before (cf. Eqn. 80 1 



and * denotes an element- by- element matrix produc 



Once again, Eqn. 90 is lower-triangular in k and requires 



only matrix multiplication and the inversion of a J-by-J 
matrix K times. G k is initialized as in Eqn. 



79 



and the 

fcth term is computed from the previous fc' < fc terms. 
The joint distribution is retrieved via the inverse trans- 
form 



jk 



(91) 



4- The \n, k n ) basis 

Expansion in either the \j, fc) or the \j, kj) basis conve- 
niently turns the master equation into a lower-triangular 
linear algebraic equation and replaces cutoffs in protein 
number with cutoffs in eigenmode number (which can 
be smaller with appropriate choices of gauge). However, 
these bases sacrifice the original tridiagonal structure of 
the master equation in the copy number of first gene, 
ri. Therefore we now consider a mixed representation, 
in which the first gene remains in protein number space 
\n), and we expand the second gene in an n-depcndcnt 
eigenbasis \k n ), 



\G) 



E 



G nk \n, k„ 



(92) 



If the rate parameter of the |fc„) basis were the regula- 
tion function q n , it would be natural to make Hq the 
entire Hamiltonian, Eqn. |64| For generality we will in- 
stead allow the rate parameter of the |fc„) basis to be an 
arbitrary n-dependent local gauge q n , such that Eqn. 63 
at steady state naturally partitions as 



= -J2 G nk \H (n) + Hi(n) 



\n,k n ) 



nk 



where 



H (n) = b+b n + pb+b m (n) , 
Hi{n) = p6+A„(n), 



(93) 



(94) 
(95) 



with b m (n) = a m - q n and A n (n) = q n - q n . Note 
that \n, k n ) is not the eigenbasis of Ho(n), but rather 
^0(^)1^7^1) retains the original tridiagonal structure in 
n, i.e. 

H (n)\n,k n ) = (g n + ri)\n, k n ) - g n \n + 1, k n ) 

-n\n— l,k n ) + pk\n, k n ), (96) 

where Eqns. [65] |67] [12] and [13] are recalled in applying 
the first term of H (n). 

Projecting the conjugate state (n, k n \ onto Eqn. 93 
yields, after some simplification (cf. Appendix C), the 
equation of motion 

5«-iG„-i,/c + ( n + l)G„ + i jfc - (g„ + n + pk)G nk 

k k 
= -ffn-l E V nfi n -i,k-i - (n+ 1) ^ V^Gn+i^-i 



+pA„G ni fe_i, 



(97) 
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where 



An 



Qn Qn i 

(-QiY 



(98) 
(99) 



and Q± = q n - 
algebraically as 



q n ±i- Eqn. 97 can be written linear 



Gh — 



(T fe ) 1 {(S-g)*diag[V^(S-G) T ] 
+(S+n) *diag[V+(S+G) T ] 
+pA*GVi}, (100) 

where X T indicates the transpose of X, * denotes an 
element-by-element product, S^ n , — 5„±i !n ' are super- 
(+) and subdiagonal (— ) matrices, T* n , = g n -i^n-i,n' + 
(n + l)<5„+i. rl ' — (g n + n + pk)5 nn > is a tridiagonal ma- 
trix, and G, which is G with the columns reversed (i.e. 
G n ( — G n ,k-t), is built incrementally in k. As with the 
\j, k) and \j, kj) bases, the fcth term is slaved to the pre- 
vious k' < k terms; in total the solution requires K in- 
versions of an N-hy-N matrix. However, here the task of 
inversion is simplified because the matrix to be inverted 
is tridiagonal. In fact, using the Thomas algorithm [40] . 
we obtain an analytic solution for the case of constant 
production of the first gene and threshold regulation of 



the second gene, as described in section pC 
The solution is initialized at k = using 



G n o 



i,(k = 0) n \G) = J2(n\n'){0n\m)p v 



Pn 



(101) 



(cf. Eqn. 26 1, where p n is known (cf. Eqn. 60 I, and the 
joint distribution is retrieved via the inverse transform 



Pnm = (n,m\G) = }^G nk (m\k„) ■ 



(102) 



Substituting Eqn. 103 into Eqn. [58] at steady state 
gives, after some simplification (cf. Appendix C), 

= 52G jk (n,m\{-H\M 

jk 

+at,g n \j,S-k n ) + a~\j,5 + k n )}, (105) 
at and a~ act as in 



where H is defined as in Eqn 
Eqns. [T2p5| and 



64 



\6±kn) = \k n±1 ) - \k n ). (106) 

We may now partition H — Ho + Hi as 

H (n) = b+b~ + pb+b-(n), (107) 
H x {n) = b+t n + pb+A n (n), (108) 

with b~ — a~ — g and T n = g — g n as in the \j, k) and 
\j, kj) bases, and 6~ (n) = a~ n — q n and A n (n) = q n ~q n as 
in the |n, k n ) basis. Noting that \j, k n ) is the eigenbasis 
of Ho(n), i.e. 



H (n)\j, k n ) = (j + pk)\j, k n ) 



(109) 



Eqn. 105 becomes, after more simplification (cf. Ap- 
pendix C) , 

= -U + pfyGj'k-^Tj-ij'Gj'k-p^Ajj'Gf^-i 
j' j' 



where Tjji is as in Eqn. |76[ 

A n' = J2ti\ n )(qn - q n )(n\f), 

n 

A tr = E0»(«+i)<" + iUXt 

n 

A w' = E0»5n-i(n-l|j')F„7, 



(110) 

(111) 
(112) 
(113) 



5. The \ j, k n ) basis 

We now consider a basis which employs both the 
constant-rate eigenfunctions \j) in the n sector and the 
rt-dependent eigenfunctions \k n ) in the m sector. Ex- 
pressing the joint distribution directly in terms of the 
eigenbasis expansion of the generating function |51j , we 
write 



and V^g is as in Eqn. 99 Linear algebraically, 

-l 



Pnm = E G jk( n , m \h k n), 
jk 



(103) 



where, as in the \j, k) and \j,kj) bases, the \j) are pa- 
rameterized by the constant rate g, and, as in the |n, k n ) 
basis, the \k n ) are parameterized by the arbitrary func- 
tion q n - The inverse of Eqn. 1 103| is 



G J k = y^,Pnm(j, k n \n,m). 



(104) 



G k 



(D + s _ r) 



k-e 



(114) 



with D k --, and S, -, defined as before (cf. Eqn. 80 1, reveal 



ing once again a lower-triangular equation (i.e. each fcth 
term is slaved to the previous k' < k terms) requiring 
only matrix multiplication and the inversion of a J-by- 
J matrix K times. Recalling Eqn. |104[ the scheme is 
initialized using 

Gj = ^2p nm (j\n)((k = 0) n \m) = ^2p n (j\n) (115) 



(cf. Eqn. 26 1 with known p n (cf. Eqn. 60 1, and the joint 



distribution is retrieved using Eqn. 103 
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B. Comparison of the representations 



The spectral representations in Sec. II A produce equa- 
tions of motion with similar levels of numerical complex- 
ity. In all cases, the the original two-dimensional mas- 
ter equation has been reduced by the lower triangular 
structure in the second gene's eigenmode number k to 
a hierarchy of evaluations of one dimensional problems. 
The bases differ in the rate parameters, or equivalently 
gauge freedoms, that one is free to choose: the \j, k) basis 
requires two constants g and q; the \j, kj) basis requires 
g and a J-valued vector qj' the \n,k n ) basis requires a 
TV- valued vector q n ; and the \j,k n ) basis requires g and 
q n - 

The bases also differ in the types of problems for which 
they are most suitable. For example, the \j, k), \j, kj), 
and |j, k n ) bases, which all expand the parent species in 
eigenfunctions \j), are best when a cutoff in j is most 
appropriate, such as when the parent distribution is a 
Poisson. The \n, k n ) basis, on the other hand, is use- 
ful when a cutoff is n is most appropriate, such as when 
the parent species is concentrated at low protein num- 
ber. Different bases are more robust to numerical errors 
for different regulation functions as well: the \n, k n ) and 
|j, k n ) bases, which both rely upon repeated manipula- 
tion of the object = q n — q n ±i, are best for smooth 
regulation functions, for which the differences between 
q n and q n +i are small; the |j, k) and \j, kj) bases on the 
other hand, which involve the deviations q—q n and q~j — q n 
respectively, are less susceptible to numerical error given 
sharp regulation functions, such as a threshold. 

As indicated in Figure [2] the \j, k) basis can be viewed 
as a special case of cither the \j, kj) basis with q~j — q (in 
which case Eqn. 87 reduces to Eqn. 75 1 or of the \j,k n ) 
basis with q n = q (in which case Eqn. 110 reduces to 



Eqn. 75 1. Although possible in principle, expanding in 
the | j, m) basis does not exploit the natural structure 
of the problem, since it neither retains the tridiagonal 
structure in n nor gains the lower triangular structure in 
k. This example explicitly shows that not all bases are 
good candidates for simplifying the master equation. 

The strength of all the spectral bases discussed in this 
section, and of the proposed spectral method in general, 
is that it allows for a fast and accurate calculation of full 
steady state probability distributions of the number of 
protein molecules in a gene regulatory network. In Fig. 
[3] we demonstrate this property for the two-gene system 
by plotting error versus computational runtime for each 
spectral basis, as well as for a stochastic simulation using 
a varying step Monte Carlo procedure |32j . For error we 
use the Jensen-Shannon divergence |41j (a measure in bits 
between two probability distributions) between the dis- 
tribution p n7n computed in the \n,m) basis (via iterative 
solution of the original master equation) and the distri- 
bution computed either via the spectral formulae in this 
section or by stochastic simulation. We plot this mea- 
sure against the runtime of each method, scaled by the 
runtime of the iterative solution in the \n,m) basis (all 
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FIG. 3: Error vs. runtime for the spectral method and 
stochastic simulation. Error is the Jensen-Shannon diver- 
gence [JT] between p nm obtained using the \n,m) basis (via 
iterative solution of the original master equation) and that ob- 



tained using the \ j, k) basis (circles; cf. Eqn. 80 1, the \ j, kj) ba- 



sis (triangles; cf. Eqn. 90 1, the |n, k n ) basis (squares; cf. Eqn. 
100 1, the \j, k n ) basis (diamonds; cf. Eqn. 114 1, or stochas- 
tic simulation [32] (dots). Runtimes are scaled by that of 
the iterative solution, 150 seconds (in MATLAB). Spectral 
basis data is obtained by varying K, the cutoff in the eigen- 
mode number k of the second gene; simulation data is ob- 
tained by varying the integration time. The input distribu- 
tion p n = n 1 e' Xl X-l/n\ + (1 - m)e~ X2 \2/n\ (from which g n 



is calculated via Eqn. 59 1 is a mixture of two Poisson distri- 



butions with Ai = 0.5, A2 = 15, and tyi = 0.5. The regulation 
function q n — g_ + (q+ — q-)n v / ' (n v + rig) is a Hill function 
with g_ = 1, q+ = 11, no = 7, and v = 2. The gauge choices 
used (cf. Fig. B are g = J2„Pngn, q = J2„Pnqn, q n = q n , 
and qj = ^2 n <j\n)q„(n\j) . The cutoffs used are J — 80 for 
the eigenmode number j of the first gene and N = 50 for 
the protein numbers n and m. Inset: The joint probability 
distribution p nm - The peak at low protein number extends 
to poo ~ 0.1. 



numerical experiments are performed in MATLAB). We 
find that the computations via the spectral bases achieve 
accuracy up to machine precision -10 3 -10 4 times faster 
than the iterative method's runtime and ~10 7 -10 8 times 
faster than the runtime necessary for the stochastic sim- 
ulation to achieve the same accuracy. Computation in 
the \ j, k) basis is most efficient since its equation of mo- 
tion is simplest (cf. Eqn. 80); the \j,kj) and \j,k n ) bases 



tend to be slightly less efficient since they require inner 
loops over £ (cf. Eqns. 90 and 1141. Fig. [3] demonstrates 
the tremendous gain in performance achieved by the joint 
analytic-numerical spectral method over traditional sim- 
ulation approaches. 
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C. An analytic solution 

In general, the equations of motion in the spectral rep- 
resentations (Eqns. 75 87 97 and 110 1 need to be evalu- 
ated numerically. In the case of the \n, k n ) basis, however, 



we can exploit the tridiagonal structure of Eqn. 97 to find 
an exact analytic solution. Specifically, in the case of a 
Poisson parent (g n = g = constant) and for threshold 
regulation, i.e. 



Qn = 



q- for n < uq 
q + for n > n a , 

setting q n — q n makes Eqn. [97] 



(116) 



gG n - 



-i,k + (n+ l)Gn+i,fc - (pk + g + n)G nk 

Hl^nnoi (H7) 



"nil! 



where 



y(-A)y 

e=i 



m .k—£i 



(118) 



(119) 



A = q + — q_, and n\ — n + I. Eqn. |117 is solved using 
the tridiagonal matrix algorithm (also called the Thomas 
algorithm [10]), as described in detail in Appendix D. The 
result is an analytic expression for the fcth column of G„fc 
in terms of its previous columns (i.e. the matrix inversion 
has been performed explicitly), 

ni(n -l)! r\ k n 
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/^«vn^(ef-i) 



1) n < no 
n > no, 



where 



N—n n-\-i 

En 



i 



- 1 

t=0 l=n ~z 
pk + g + n n* 
gn 

n . 

\ -» n! 



Vn-l 



7-9 



1=0 



^— ' i!(n — ill" 



(120) 

(121) 
(122) 
(123) 
(124) 
(125) 



and iV is the cutoff in protein number n. Along with the 
analytic form of the mixed product 



(ml fc, 



(-\) k e-^ql l k\ 

min(m,fc) 

x E 



£=0 



£\(m-£)\(k-ey.(-q n y 



(126) 



(cf. Appendix A), Eqn. 120 in the limit N — > oo consti 



tutes an exact analytic solution for the joint distribution 
p„ m , as calculated using Eqn. |102| 



D. The threshold-regulated gene approximates the 
on/off gene 

If a gene is regulated via a threshold function (cf. Eqn. 



116), its steady state protein distribution p m can be well 



approximated by the two-state process discussed in Sec. 
|I C| To make the connection clear, we first observe that 
the off-state (z = — ) corresponds to the first gene ex- 
pressing the same or fewer proteins n than the threshold 
no, i.e. 



Pn 



= E p« 



(127) 



and the on-state (z = +) corresponds to the first gene 
expressing more proteins than the threshold, i.e. 



Pn 



= E P« 

n>n a 



(128) 



The dynamics of p^ are then obtained by summing the 
master equation for two-gene regulation, Eqn. |58| over 
either all n < no or all n > no, giving 

Pm = P [l±Pm-l + ( m + l)Pm+l ~ (<7± + m )Pm] 

Tnip nim ± (129) 
where n% = no + 1. Making the approximations 



Pm 
17- 

7T+ 



p(m\-) 
p(m\+) 



p(m\n ) 
p(m\ni) 



Pn 



Pn 
Pn\m 



where 



= E p ™ = E p "' 

m n<no 

= E p ™ = E p » 

m n>riQ 



(130) 
(131) 

(132) 
(133) 



are the total probabilities of being in the off- and on- 
states respectively, and noting from Eqn. [59] that g no = 



nip ni /p no , Eqn. 129 at steady state becomes 

= QzPm-i + i m + l)Pm+i - {iz + m)p z m 

+ ^n zz ,p z m , (134) 



with z — ± and 



n = 



(135) 



12 



where 



nip ni 



(136) 



Eqns. p4p35] have the same form as Eqns . [Tj and |47| at 
steady state with n — ► m and o — > 5, and Eqn. |l36| relates 
the effective switching rates w± to input and regulation 
parameters p ni , tt±, and ni, and the ratio p of the degra- 
dation rate of the second gene to that of the first. Note 
that Eqn. 136 satisfies 



7T_ 
7T+ 



(137) 



in agreement with Eqn. |48[ and exhibits the intuitive 
behavior that increasing p (i.e. decreasing the relative 
response rate of the first gene) is equivalent to decreasing 
the switching rates u>±. 

A comparison of the distributions of a threshold- 
regulated gene with those of an on/off gene for vari- 
ous parameter settings reveals that Eqns. |130|131| are 
a good approximation. Fig. [4] shows a demonstration for 
a threshold-regulated system with a Poisson input distri- 
bution. In the first column, the mean g of the input lies 
above the threshold hq , making the output more likely to 
be in the on-state, i.e. 7r + > 7r_; in the second column, 
g < rio, making 7r + < 7r_. In the first row p < 1; in 
the second row p > 1, corresponding to lower effective 
switching rates lo± and producing bimodal distributions 
with peaks near the on/off rates q±. In all examples, 
the approximation as a two-state process with switch- 
ing rates given by Eqn. |136| agrees well with the actual 
output from threshold regulation. 



III. REGULATION WITH BURSTS 

The final system we consider combines the multi-state 
process used to model bursts of expression in Sec. [TJwith 
gene regulation as discussed in Sec. |TT] Specifically we 
consider a system of two species, with protein numbers 
n and m, existing in Z possible states, distinguished by 
the settings of the two production rates g z and q z re- 
spectively, where 1 < z < Z. Regulation is achieved by 
allowing the rates of transition among states affecting the 
production of the second gene to depend on the number 
n of proteins expressed by the first gene. Recalling Eqns. 
[T] and [58] the master equation describing the evolution of 
the joint probability distribution p^ m reads 

Pnm = 9zPn-l, m + (™ + l)Pn+l,m ~ 0* + ™)Pnm 

+P [izP^m-i + (m+ l)K,m+i - (iz + m)p z nm ] 

+E fi »'(^' ( 138 ) 

z' 

where the dependence of the stochastic matrix Q, zz > on n 
incorporates the regulation. 

As with the previously discussed models, Eqn. |138| 
benefits from spectral expansion, and for simplicity we 



p = 0.1 




o 10 20 30 
number of proteins, m 



10 20 30 
number of proteins, m 



FIG. 4: Protein distributions for a gene regulated by a thresh- 
old function (dots; calculated via Eqn. 75 1 and a gene with 
two stochastic states (circles; calculated via Eqn. 43 1. The 
relationship between regulation parameters and state transi- 
tion rates is given by Eqn. |136| In all panels the input to 
the regulation is a Poisson distribution with mean g = 7, 
and the regulation rates (cf. Eqn. 116 1 are g_ = 2 and 
q+ = 15. In the first column the threshold is no = 4 making 
7r+ = 0.827 > 7T„ = 0.173; in the second column no — 8 mak- 
ing 7r+ = 0.271 < tt- — 0.729. In the first row the ratio of the 
second gene's degradation rate to that of the first is p = 0.1; 
in the second row p — 10. 



present only the formulation in the \j, k) basis, parame- 
terized by constant rates g and q respectively, as in Sees. 
|IB| and |II A2| As before the first step is to define the 
generating function 



\G Z 



Pnm\ n , m ), 



(139) 



with which Eqn. |138[ upon summing over n and m 
against |n, m), becomes 

\G z ) = -H z \G z ) + J2&zz'\G z ,), (140) 

z' 

where 

(141) 
(142) 
(143) 
(144) 
(145) 

and £l zz > is £l zz '(n) with every instance of n replaced by 
the number operator a^a~. Defining b~ = a~ — g and 
&~ = — q, we partition the Hamiltonian as H z = 
H + Hf, with 



H z 


= b+b- z + pb+b 




= a+ - 1, 




= a m - !. 


Kz 


= K ~ 9z, 




= a m - q z , 



Ho = h tb n + pb+b m 



(146) 
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the operator of which \j, k) is the eigenbasis, i.e. 

H \j,k) = (j + pk)\j,k), (147) 

and 

Hi = b+T z + pb+A z (148) 

capturing the deviations T z = g — g z and A 2 = q — q z of 
the constant rates from the state-dependent rates. Upon 
expanding the generating function in the eigenbasis, 



\Gz) = Y / G z jk \j,k), 

jk 



(149) 



and taking dummy indices j — > f and k — > k' , projecting 
the conjugate state (j, k\ onto Eqn. 140 gives 



Gjk — —(j + pk)G z ]k — F z Gj_ lk — A z Gj k _ x 



(150) 



where the components of Cl zz i need only be evaluated in 
the j sector, not the k sector, because the transition rates 
depend on only n, not m (cf. Eqn. 1381. Like Eqns. 43 
and |75| Eqn. |150| is subdiagonal in k and thus far more 
efficient to solve than the original master equation, Eqn. 
138[ as we demonstrate for a special case in the next 
section. The joint distribution is retrieved from G| fc via 
inverse transform, 



Pnm = J2(n\j)Gj k (m\k), 
jk 



(151) 



with the mixed products calculated as in Appendix A. 



A. The four-state process 



As a simple example of the model in Eqn. |138| we con- 
sider a system in which each of the two species has an 
on-state and an off-state, and the transition rate of the 
second species to its on-state is a function of the number 
of copies of the first species. This system models both 
(i) a single gene for which the production of proteins de- 
pends on the number of transcripts, and each is produced 
in on- and off-states by the binding and unbinding of ri- 
bosomes and RNA polymerase respectively, and (ii) one 
gene regulating another with each undergoing burst-like 
expression. 

There are a total of Z = 4 states, i.e. 



P'n 



(Pnm ' Pnm ' Pnm, ' Pnm ) ' 



(152) 



where the first signed index denotes the state of the first 
gene (with protein count n) and the second signed in- 
dex denotes the state of the second gene (with protein 
count m). Defining g± as the production rates of the 



first species in its on- (+) and off-states (— ), and simi- 
larly q± for the second species, the production rates of 
the Z = 4 states are: 



9z = {9-,9+,9-,9+), 
Qz = {q-,q-,q+,q+). 



(153) 
(154) 



Defining u± as the transition rates of the first species to 
(+) and from (— ) its on-state, and similarly a± for the 
second species, the transition matrix takes the form 



uj + -a + (n) 


u>_ 









—lu- — a+{n) 





OL- 







— lu + — a- 


UJ- 





a+{n) 




— LU- — 



(155) 

The simple form a + (n) = cn v for constant c and integer 
v corresponds to the first species activating the second 
as a multimer, with v the order of the multimerization. 
In the limit of fast switching this description reduces to 
a Hill function with cooperativity v [28 . Recalling that 
t>+ = a+ — 1 and b~ = a~ — g, the n-dependent terms of 
U\&zz' li') are evaluated as 

(j\a+(atK)\f) = c(j\[(Pt + l)(5n + 9)]"\f), (156) 

Since 6+ and b~ raise and lower \f) states respec- 
tively (cf. Eqns. 22p3 l, the modified transition matrix 
(j\fozz'\j') is nearly diagonal, with nonzero terms only 
for \j -fj < v. 

Eqn. |150| at steady state, 



-U + P k)G% - r z G*_ lifc + ^zZ(j\^'\f) G h 

z' j' 

= A,q, fe -i, (157) 

is solved successively in k, requiring the inversion of a 
4J-by-4J matrix K times. It is initialized at k = 
by computing the null space of the left hand side and 
normalizing with ^2^00 = 1 (°f- Eqn. 35 1. The joint 
distribution pf lm is retrieved via inverse transform (Eqn. 
15ll). 



With v = 2, a typical solution of Eqn. |157| takes a few 
seconds (in MATLAB), which, depending on the cutoff 
N, is ~10 2 -10 3 times faster than direct solution of the 
master equation, Eqn. 138 by iteration, for equivalent 
accuracy. The advantage of such a large efficiency gain 
is that it allows repeated evaluations of the governing 
equation, necessary for parameter inference or optimiza- 
tion [34]. We demonstrate this possibility in the next 
section by finding and interpreting the solutions that op- 
timize the information flow from the first to the second 
species. 
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B. The information-optimal solution 

Cells use regulatory processes to transmit relevant in- 
formation from one species to the next j42l H3J 04j |45j 06] . 
Information processing is quantified by the mutual infor- 
mation /, which, between the first and second species in 
the four-state process, is 



Pnm 
PnPn 



(158) 



where the distributions p„ m , p n , and p m are obtained 
from summing the joint distribution p nm (cf. Eqn. 151 1, 
and the log is taken with base 2 to give / in bits. 

Upon optimization of I for the four-state process, two 
distinct types of optimal solutions become clear: those 
in which the distribution p nm has one peak, and those in 
which p nm has two peaks. The former occur when copy 
number is constrained to be low, and switching rates are 
constrained to be near the decay rates of both species, 
producing a single peak at low copy number (see lower 
left inset of Fig. |5|3). As these constraints are lifted, it is 
optimal for the switching rates of the parent species to 
become much less than the decay rate. The slow switch- 
ing produces a second peak whose location is specified by 
the on-rate of each species (see upper right inset of Fig. 

§3). 

To quantify the transition between the two types of 
solutions, we numerically optimized mutual information 
over parameters g+ 1 q+, lu~, uj + , a_ , and c (the off-rates 
g_ and q_ were fixed at 0; the cooperativity v was fixed 
at 2; and the decay rate ratio p was fixed at 1). In- 
formation may always be trivially optimized by allowing 
infinite copy number or arbitrary separation of relevant 
timcscalcs. We limit copy number by constraining the 
gain 



7 



r + A 



(159) 



defined as the average of the parent gain T = g + — g_ 
and the child gain A = q + — q_. Since g_ = g_ = 
0, the maximum number of particles is dictated by the 
on-rates g + and q + , and thus constraining 7 limits the 
copy number. We limit separation between the switching 
timcscalcs and the decay timescales by constraining the 
stiffness 



1 



|log 10 w_| + |log 10 w + | 



+ |log 10 a_| + |lo gl0 [a + <n")]|), (160) 

where the average (n v ) is taken over p n . Stiffness a is the 
average of the absolute deviation of (the logs of) all four 
switching rates from the (unit) decay rates, so constrain- 
ing cr prevents fast or slow switching. Gain is fixed by 
varying g + and q+ such that 7 is a constant, and stiffness 
is constrained by optimizing the objective function 



I - Act 



(161) 



for a given value of the Lagrange multiplier A. 

As shown in Fig. |5|V, one-peaked solutions are more 
informative at low stiffness, while two-peaked solutions 
are more informative at high stiffness. We compute the 
convex hulls of the one- and two-peaked data to remove 
suboptimal solutions, and the transition occurs at the 
stiffness value at which the convex hulls intersect (cf. 
Fig. [5|\). Repeating this procedure for many choices of 
gain allows one to trace out the phase transition shown 
Fig. [5)3, which makes clear that one-peaked solutions are 
most informative at low stiffness, two-peaked solutions 
are most informative at high stiffness, and the critical 
stiffness decreases weakly with increasing gain. 



IV. CONCLUSIONS 

The presented spectral method exploits the linearity of 
the master equation to solve for probability distributions 
directly by expanding in the natural eigenfunctions the 
linear operator. We demonstrate the method on three 
models of gene expression: a single gene with multiple 
expression states, a gene regulatory cascade, and a model 
that combines multi-state expression with explicit regu- 
lation through binding of transcription factor proteins. 

The spectral method permits huge computational 
gains over simulation. As demonstrated for all spec- 
tral expansions of the two-gene cascade (cf. Fig. [3]), di- 
rectly solving for the distribution via the spectral method 
is ~10 7 -10 8 times faster than building the distribution 
from samples using a simulation technique. This massive 
speedup makes possible optimization and inference prob- 
lems requiring full probability distributions that were not 
computationally feasible previously. For example, by op- 
timizing information flow in a two-gene cascade in which 
both parent and child undergo two-state production, we 
reveal a transition from a one-peaked to a two-peaked 
joint probability distribution when constraints on protein 
number and timescale separation are relaxed. We empha- 
size that this optimization would not have been possible 
without the novel efficiency of the spectral method. 

The spectral method also makes explicit the linear al- 
gebraic structure underlying the master equation. In 
many cases, such as in two-state bursting and the two- 
gene threshold regulation problem, this leads to analytic 
solutions. In general, such as shown in the case of the lin- 
ear cascade, this leads to a set of natural bases for expan- 
sion of the generating function and reveals the features 
of each basis that are better suited to different types of 
problems. Specifically, bases in which the parent species 
is expanded in eigenfunctions are best when the parent 
distribution is Poissonian, and bases in which the parent 
is left in protein number space are best when the par- 
ent distribution is concentrated at low protein number. 
As well, bases in which the eigenfunctions of the child 
depend on the number of copies of the parent's protein 
are best suited for smooth regulation functions, whereas 
a basis in which the eigenfunctions of the child are pa- 
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0.4 0.6 0.8 

stiffness, a 



FIG. 



5: A: Mutual information / (cf. Eqn 
ness a (cf. Eqn. |160[) for fix ed ga in (7 = 



161 



158 1 versus stiff- 
16, cf. Eqn. [159} , 
for A values between 10 -3 



obtained by optimizing Eqn. 

and 10 1 . Squares denote solutions whose joint distribution 
p nrH has one peak (cf. B, lower left inset), and dots denote so- 
lutions for which p nm has two peaks (cf. B, upper right inset). 
Solid lines show the convex hulls of the one- and two-peaked 
solutions. Dotted lines indicate the stiffness value at which 
the hulls intersect and the stiffness values of the hull points 
to the left and right of the intersection. B: Phase diagram 
between one- and two-peaked optimal solutions in the gain- 
stiffness plane. Circles and left and right error bars at each 
gain are determined by the stiffness values at the intersection 
of the one- and two-peak convex hulls and at the hull points to 
the left and right of the intersection respectively (see dotted 
lines for the example case in A) . Solid line shows a line of best 
fit. Insets show examples of one- (lower left) and two-peaked 
(upper right) optimal distributions p„ m . 



rameterized by a constant is more numerically robust for 
sharp regulation functions such as thresholds. In all cases 
the linear algebraic structure of the spectral decomposi- 
tion yields numerical prescriptions that greatly outper- 
form simulation techniques. We anticipate that the com- 
putational speedup of the method, as well as the removal 
of the statistical obstacle of density estimation inherently 



limiting simulation-based approaches, will make spectral 
methods such as those demostrated here useful in ad- 
dressing a wide variety of biological quesitons regard- 
ing accurate and efficient modeling of noisy information 
transmission in biological systems. 



APPENDIX A 

In this appendix we describe two ways to compute the 
mixed products (n\j) and (j\n) between the protein num- 
ber states |n) and the eigenstates by direct evaluation 
and by recursive updating. 

The direct evaluation follows from Eqns. [9j [8] [20] and 
[TO] and the fact that repeated derivatives of a product 
follow a binomial expansion. Introducing g as the rate 
parameter for the \j) states, 



Hi) - 



dx 



n\x){x\j) 



2m 

dx e^ x - l \x-l)i 
2m x n+1 



n 
1 



n! ^ ii(n-ey. 



x=0 



x=0 



e=o v 



r - iy_ 
1 



x=0 

[9 n ~ 



J 



(-l) j e-9g n M nj , 



where 



min(n,j) 



^ t\(n-iy.(j-£)\(- q y 



Similarly, noting Eqns . [21] and [7] 

(j\n) = n\(- g y£ nj , 



(162) 

(163) 
(164) 

(165) 

(166) 
(167) 

(168) 

(169) 



with £ n j as in Eqn. 168 Eqns. 169| and 167 clearly reduce 
to Eqns. |26p7| for the special case j = 0. 

It is more computationally efficient to take advantage 
of the selection rules in Eqns. T2][T5 and 22"]|25" to compute 
the mixed products recursively. For example, using Eqns. 
[221 [161 andflil 



(n\j + l) - (n\b+\j) = (n\(a+-l)\j) 
= (n - - (n\j), 



(170) 



which can be initialized using (n|0) = e 9 ' g n ln\ (cf. Eqn. 
167) and updated recursively in j. Eqn. 170 makes clear 



that in n space the (j + l)th mode is simply the (negative 
of the) discrete derivative of the jth mode. 
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Alternatively, Eqns. 15 T7] and 23 give 



(n + l)(n+l\j) = {n\cT\j) = {n\(b~ + g)\j) 

= 3{n\j-l) + g{n\j), (171) 

which can be initialized using (0|j) = (— l) 3 e~ 9 (cf. Eqn. 
1671 and updated recursively in n. 



One may similarly derive recursion relations for (j\n) 



(j\n+l) = (j-l\n) + (j\n), (172) 
(j + l)(j+l|n) - n(j\n-l)-g(j\n), (173) 

initialized with (j|0) = (—g) 3 /j\ or (0|n) = 1 respec- 
tively (cf. Eqn. 169| and updated recursively in n or j 
respectively. 

One may also use the full birth-death operator b + b 
to derive the recursion relations 



(n + l)(n+l|j) - (g + n-j)(n\j)-g(n-l\j), (174) 
g(j\n+ 1) = (g + n- j)(j\n) -n(j\n- 1), (175) 



initialized with (0|j) 



(-l) J e" 9 and = 



(-l) 3 e-"{g-j) (cf. Eqn.|l67[), and (jjO) = {-g) 3 /j\ and 
til 1 ) = (-5) J (1 ~ j/g)/WT c i- Eqn. [169]), respectively, 
and updated recursively in n. We find Eqns. 174|175| are 



more numerically stable than Eqns. |170||173 



as the for- 
mer are two-term recursion relations while the latter are 
one-term recursion relations. 



or, noting Eqn. 55 and the fact that T(z + 1) = zT(z) 
for any z, 



G{x) = J2 



^+ r 0' + ^+ + 1 ) 
lo+ + Ld- r(cj + + 1) 



E 



W_ T(j + UJ + ) 

r(i + w+ + cj_ + i) j!' 

w + (i + ^+)r(j + w+) 



E 



cj + + cj_ r( w +) 

(w + + w_)r(w + + w_) 
0' + w + + w_)r(j + w + + w_) j! 
r 0' + ^+) r(^ + + ^_) ^ 
r(w+) r{j + uj + + uj-) j\ 



= $[uj + ,uj+ + uj_;y], 

as in Eqn. [56] 

The marginal p n is given by 



(179) 



(180) 

(181) 
(182) 



(183) 



(cf. Eqn. 10). Using Eqn. 56 and the derivative of the 
confluent hypergeometric function, 

d^[a, /?; y] = J} P \, <f[a + n,/3 + n;y], 



T(a) T(n + [3) 



one obtains Eqn. [57] 



(184) 



APPENDIX B 



In the limit <?_ = 0, Eqn. 54 reads 



G(x) = ^e"$[u;_,w + +w_ + l;-y] 

-$[w + ,w + +w_ + l;y] > (176) 



LU + + LU- 



where ?/ = e a +^ x ^ . Using the fact that [17 
e^[a ! /3;-t/] = $[/3-a ! /3;j/] ! 
Eqn. [176] can be written 



(177) 



G(x) 



-$[uj + + + l;y] 

-$[w + ,w++w_ + l;y], (178) 



APPENDIX C 



In this appendix, we fill in the details of the derivations 
of the equations of motion for the latter three of the four 
spectral bases discussed in Sec. |II A| 



The | j, kj) basis 



Projecting the conjugate state (j, kj\ onto Eqn. 83 (in 
which dummy indices j and k are changed to f and k' 
respectively) gives 



= ]>>■' + pfcW)(fc#>>G 



j'k' 



j'k' 



■J2(j\b^n\f)(kjWr)G. 



j'k' 



j'k' 



+pY.( k M\kf)(j\&n(j')\j')G 3 ' k ,. (185) 

j'k' 
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From the orthonormality of states, the first term of Eqn. 
|185| simplifies to 

J2U + phf){kj\^G jV = (j + pk)G jk . (186) 



Recalling Eqn. 30 the product {kj\k'y) simplifies to 

\k—k' 



(k 3 \k'y 



(k-k')\ 



-6(k-k' + l), 



(187) 



with Qjji = qj — qji, whereupon Eqn. 185 separating the 
part of its second term which is diagonal in k from that 
which is subdiagonal and applying Eqn. [24] to its third 
term, becomes 

= (j + pk)G jk 



3' 

k'<k j' 
k'<k j' 



i-Qj/ >l 1 

(k-k')l 



Gj'k> 



33' ' 



{ -Q33') k - k '- 1 

(k-k'~ 1) 



Gyv, (188) 



, 3 



with Tjji as in Eqn. 76 and 



A lf = (j\A n (j')\j') = (j\(q f - q n )\f) (189) 



ti\(qj-qn)\f) 
E^'l n )fe' ~ 1n){n\f) 



(190) 
(191) 



(where the orthonormality of \j) states is used in going 
from Eqn. 189 to Eqn. 190 1. Defining I = k — k' and 



yt _ ( Qjj') 

33' l\ 



(192) 



Eqn. |188| can be written more compactly as Eqn. [87] 



The \n, k n ) basis 



Projecting the conjugate state (n, k n \ onto Eqn. 93 (in 
which dummy indices n and k are changed to n' and k' 
respectively) gives 

= E(-9«' +n' + pk')(n\n')(k n \k' nl )G n >k' 

n'k' 

-^2g n ,(n\n' + l){k n \k' n ,)G n , k , 

n'k' 

-Y,n'(n\n' -l)(k n \k' n ,)G n , k , 

n'k' 

+ P J2&n'(n\n')((k-l) n \k' n ,)G n , k >, (193) 



where A n — q n — q n . Noting that, as in Eqn. 30 



(k n \K l±1 )= ( { ^_ kf)l e(k-k' + i), 



(194) 



where = q n — q n ±i, Eqn. 193 becomes 
= (g n + n + pk)G nk 

— TWt ^n — l.k' 



k'<k 



(k - k')\ 



(_rt+\k-k' 

k'<k v ' 
-pA n G nik -i. 



(195) 



Separating the parts of the second and third term that 
are diagonal in k and defining i = k — k' and 



(~Qn) 1 



(196) 



Eqn. [195] becomes Eqn. [97] 



The | j, k n ) basis 



Substituting Eqn. |103| into Eqn. [58] at steady state 
gives 

= G jk {g n ^i{n ~ l,m\j, k n -i) 
jk 

+(n + l)(n+ l,m\j, k n+1 ) - (g n + ri)(n, m\j, k n ) 
+p [q n (n,m- l\j,k n ) + (m + l)(n, m + l\j, k n ) 
-(q n + rn)(n, m\j, k n )]} , (197) 

or, in terms of raising and lowering operators (cf. Eqns. 
T4fl5] , 

= y^Gjfc(Tt,m| {atg n \h k n -i) 
jk 

+a~\j,k n+ i) - (g n + a+a~)\j, k n ) 
+P [ati9n\j, 

Kfij ~r & m \]-} &n) 

-(g„ + a+a-)|j,fc„)]}. (198) 
Using the definitions in Eqns. [64|68] Eqn. [198] can be 



written as Eqn. 105 

Using Eqns. |107fll09] Eqn. |105| can be written 



= - 



hm\Y,(f + pk')\f,K)G 



j'k> 



j'k< 



i,m\^2b+t n \f,k' n )G fk/ 



j' k < 



i\pJ2^n-qn)(n\f)b+\k' n )G f 



j'k< 



Y.gn-xin-lW^-k'^G^ 



ji k > 



Y,{n+l)(n + l\j')\8 + k' n )G rk ,, (199) 
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where 



and the solution is obtained by backwards iteration with 



\S±k n ) = \kn±i) ~ \k„), 



(200) 



and the dummy indices j and k have been changed to j' 
and k' respectively. Using Eqn. 194 to note that 



(k n \S±k' n 



{ -^ r 0(k-k % (201) 



where Q„ = q n — q n ±i, we multiply Eqn. 199 by (k 



and sum over m to obtain 

= -(n\J2U' + pk)\f)G 3 > k 
j' 

j' 



-P^2iQn - qn){n\f)Gj',k-l 
j' 

k 

-^9n-i<n-l|/)^y-^, fe _, 



^(n+mn+llf^V+Gj,^, (202) 



e=i 



in which we exploit the completeness of 
E ra \m)( m\ = 1, and is as in Eqn. 



m) states, i.e. 
Multiply- 



On 



ing Eqn. 202 by (j\n), summing over n, and exploiting 
l n )( n l = 1 for the first two terms, we obtain Eqn. 

nni 



APPENDIX D 

In this appendix we explicitly solve for G nk in Eqn. 

, algo- 



117 using the tridiagonal matrix, or Thomas 



rithm. We start by identifying the subdiagonal, diagonal, 
superdiagonal, and right hand side elements of Eqn. |117| 
respectively, as 



G h N = R' N , (211) 
GLi = K-i - C' n -iG k n (n = N . . . 1) (212) 

(where k has been moved from subscript to superscript 
for ease of reading) . 

Computing the first few terms of Eqn. |208| reveals the 
pattern 



C' n = -(n+l)- k . 

Vn+l 



(213) 



where 



i-l 



ti = T,min y9 n - l E(p k +^ ( 214 ) 

i=0 ' v 1=0 

with the convention that J\ b a [-] = 1 if a > b. Note that 
since Y^i= Q {pk+l) = F(pk+i)/T(pk), we may also use the 
integral representation of the Gamma function to write 

re! 



u Ul ^T ff""' f dte-H^-V (215) 
±-r C dte-H^-^Y- 



_ -g n ~H* (216) 

-z!(n-z)! 



= fW)Jo dte ~ tt(Pk ~ 1){9 + t)r 
Using Eqn. |210| it is immediately clear that 



R"n<n ~ 0- 



A n = g 

B n = ~{pk + g-\ 

C n =71+1 



(n = l...N), (203) 

(n = 0...N), (204) 

(re = 0...N- 1), (205) 

(n = 0...N), (206) 



where N is the cutoff in protein count re and rei = reo + 1. 
Auxiliary variables are defined iteratively as 



The first nonzero term is 

R > = nifik Vn 1 

"° arm n k f k — 1 ' 
y ,L '/n -l "o 1 

where we have defined 

k = pk + g + n rfc 
9^ Vn-l 

Further iteration of Eqn. |210| makes clear that 

'l V? 1 \ 



(217) 



(218) 



(219) 



(220) 



C' 



i?n = 



Co 
So' 



R, 



"n>no 



B n C n _-yA< n 

-Rn R"n— \Af 



B,, 



G' n -\An 



(n= 1 



(n = 1 



(207) 
JV- 1), (208) 

(209) 
.N), (210) 



»i 



^ t^no \ nil 



i K - 1)! 

no— ^ i=n 



g re! ?7™ 



(221) 

1 

-All^ ( 222 ) 

1 ;-«„ e i 1 



where 



/*=^-^^(4-l)4V (223) 

"1 '/ra 
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Computing the first few terms of Eqn. 212 reveals that 

(224) 



(~ik I k i \ rtl pfe 

>-^„^„„ — \t n L ) Jx n>n r n 



n>n 



m (n - 1)! rfc 
9 n\ n k nn _ 



Jk i 

Vn f pt rr 1 

k Jk r n J_J_ fc _ 1 

rao-1 i =r ,„ i 



(225) 



These results are summarized in Eqns. 120P25 



where 



i=0 l=n ^ 



At the threshold Eqn. 212 gives 



no 



1 



(226) 



(227) 



and since R' n<n<) — 0, the solution is easily completed 
using Eqn. |212| giving 



i=l 

9 n\ j)* 



f k _ 1 

n -l "0 



(228) 



(229) 
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